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ABSTRACT 

Dust infrared emission possesses scaling properties that yield powerful results 
with far reaching observational consequences. Scaling was first noticed by Rowan- 
Robinson for spherical shells and is shown here to be a general property of dust 
emission in arbitrary geometries. Overall luminosity is never an input parameter of 
the radiative transfer problem, spectral shape is the only relevant property of the 
heating radiation when the inner boundary of the dusty region is controlled by dust 
sublimation. Similarly, the absolute scales of densities and distances are irrelevant; 
the geometry enters only through angles, relative thicknesses and aspect ratios, and 
the actual magnitudes of densities and distances enter only through one independent 
parameter, the overall optical depth. That is, as long as the overall optical depth 
stays the same, the system dimensions can be scaled up or down by an arbitrary 
factor without any effect on the radiative transfer problem. Dust properties enter only 
through dimensionless, normalized distributions that describe the spatial variation of 
density and the wavelength dependence of scattering and absorption efficiencies. 

Scaling enables a systematic approach to modeling and classification of IR spectra. 
We develop a new, fully scale-free method for solving radiative transfer, present exact 
numerical results, and derive approximate analytical solutions for spherical geometry, 
covering the entire range of parameter space relevant to observations. For a given type 
of grains, the spectral energy distribution (SED) is primarily controlled by the profile 
of the spatial dust distribution and the optical depth — each density profile produces 
a family of solutions, with position within the family determined by optical depth. 
From the model SEDs presented here, the density distribution and optical depth can 
be observationally determined for various sources. 

Scaling implies tight correlations among the SEDs of various members of the 
same class of sources such as young stellar objects, late-type stars, etc. In particular, 
all members of the same class occupy common, well defined regions in color-color 
diagrams. The observational data corroborate the existence of these correlations. 

Key words: infrared - dust ~ radiative transfer - stars: late-type - young stellar 
objects 



1 INTRODUCTION 

For many astronomical objects, the radiation we receive has 
undergone significant processing by surrounding dust, and 
its detailed interpretation requires considerable theoretical 
effort. Input to the necessary radiative transfer calculations 
includes properties of the original radiation, optical charac- 
teristics of individual dust grains, and the dust density dis- 
tribution. Specifying all the relevant properties traditionally 
involved a rather large number of input parameters, creat- 
ing two major practical problems. First, the volume of pa- 



rameter space that must be searched to fit a given set of 
observations can become prohibitively large. Second, and 
more serious, even when a successful fit is accomplished, its 
uniqueness is questionable and the model parameters can- 
not be trusted as a reliable indication of the source actual 
properties. 

In fact, much of this input is redundant because the 
radiative transfer problem possesses general scaling proper- 
ties that drastically reduce the number of independent input 
parameters. Rowan-Robinson (1980, RR hereafter) was the 
first to utilize scaling in his extensive study of IR emission 
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from spherical dust shells. In his formulation of the radia- 
tive transfer problem, the overall luminosity, as well as other 
quantities, never enters. The full power of this scaling was 
exploited in our previous work on late-type stars (Ivezic & 
Elitzur 1995, IE95 hereafter). There we show that, given the 
dust optical properties, for each circumstellar shell both the 
dynamics and IR emission are successfully described by a 
single parameter - the overall optical depth. This explains 
the many correlations displayed by the data for late-type 
stars. Here we derive the general form of these results, ex- 
tending the scaling analysis to arbitrary geometries and den- 
sity distributions. Formal discussion of the general scaling 
properties of the radiative transfer equation is presented in 
Section 2, and applied in Section 3 to radiatively heated 
dust in arbitrary geometries. In Section 4, more detailed, 
concrete results are derived for spherically symmetric sys- 
tems, leading to a new, exact method for solution of the 
radiative transfer problem that can be readily extended to 
arbitrary geometries. Thanks to scaling, the effects of the 
independent input properties on emerging radiation can be 
searched systematically. In Section 5 we present such an ex- 
haustive search for spherical systems and the entire relevant 
range of input properties. We explore the effects of various 
density distributions, grain size and chemical composition, 
sublimation temperature, and spectral shape of the external 
radiation. Section 6 contains a summary of our results and a 
discussion of observational implications. A subsequent paper 
will present detailed analysis of the IRAS data base, which 
displays strong evidence for the scaling properties derived 
here. 



2 GENERAL SCALING PROPERTIES 

The equation of radiative transfer in steady-state is 



— iix{Sx — Ix 



(1) 



dix 

M 

Here kx is the overall extinction coefficient at wavelength A 
including both absorption and scattering, ka = KaA + «sA, 
and the corresponding optical depth element along the path 
is drx = nxdi. We will also make use of the albedo uix = 
«:sA /ka . In addition to boundary conditions, solution of this 
equation for the intensity Ix requires two input functions at 
every point r - the source function Sx (r) and the absorption 
coefficient «:A(r). The equation involves quantities with di- 
mensions of either intensity or length. The intensity Ix and 
source function Sx belong to the first type, the path-length 
£ and absorption coefficient kx, which has dimensions of in- 
verse length, to the second. Therefore the formulation of any 
radiative transfer problem can involve only two dimensional 
scales, one for intensities the other for lengths. All other rel- 
evant scales must be dimensionless, involving ratios with the 
two dimensional scales. 



2.1 Length Scales 

Consider first the length-type quantities. Because all dis- 
tances can be replaced by equivalent dimensionless ta, the 
radiative transfer problem is independent of linear size 
scales. The only geometrical properties that can enter the 



problem involve angular sizes and the spatial variation of op- 
tical depth. To further explore this scaling, introduce some 
arbitrary length scale ri. Any position r = {r,0,(j)) can be 
specified instead by the dimensionless vector y = {y,8,4)), 
where 



y = r/ri 



(2) 



is dimensionless distance from the (arbitrary) origin. De- 
note a radiative transfer path by V and positions along this 
path by (y\'P), where yr\ is distance from the path's clos- 
est approach to the origin. Then the path-length element is 
dl = r\dy and the radiative transfer equation (0) becomes 



dIx{y;T) 



rl{V)n{y;V)[Sx{y;V) ~ Ix{y;V)]. 



(3) 



dy 

Here rJ^iP) is the total optical depth along the path at 
wavelength A and 

1 drx{y;V) i^x{y\V) 



rI{V) dy J tixiy;V)dy 



(4) 



is a dimensionless function that describes the spatial vari- 
ation of opacity. In general 77(3/; V) depends also on wave- 
length. This dependence disappears when the properties of 
individual absorbers do not vary with position, in which 
case ri{y) = ni,{y)/J ns,(y)dy, the normalized dimension- 
less density distribution of absorbers. If in addition the ab- 
sorber abundance does not vary, i.e., ria/n is constant, rj 
becomes simply the normalized dimensionless density dis- 
tribution n{y)/J n{y)dy. 

In its new form, the radiative transfer equation (^ no 
longer contains any quantities with dimensions of length; 
physical dimensions enter only as ratios, i.e., angles and as- 
pect ratios. Furthermore, the length scale ri has disappeared 
altogether, thus it is not required for a complete solution. 
This quantity enters into the problem only through the prod- 
uct riKx = Tx rj. Therefore, any scale associated with it is 
contained in the overall optical depth because the func- 
tion rj sets its own scale from the normalization 



v{y;'P)dy = i 



(5) 



for any path V. This shows that radiative transfer is scale 
invariant. The physical dimensions of any system can be 
scaled up and down arbitrarily without any effect on its 
radiative properties as long as optical depths and relative 
matter distribution remain the same. Two systems with en- 
tirely different dimensions and absorption coefficients but 
the same overall optical depths will produce the same radia- 
tive intensities if they have self-similar distributions of opac- 
ities and source functions. Since the length scale ri does not 
enter into the radiative transfer problem, a corollary is that 
it can never be determined from the solution. The only way 
to determine actual length scales is through some indepen- 
dent determination of the size of the system or the distance 
to it. 

The scale-invariance of radiative transfer refiects the 
fact that, in principle, the problem could be formulated 
using optical depth as the independent variable, differen- 
tiating with respect to dTx{y;'P) = tJ^ {'P)rj{y;V)dy. In 
that case, location in space (the argument of Ix and Sx) 
is specified by its optical depth along the chosen path, 
Tx{y;P) ~ r^{P) J^ri{u;P)du. This scale-invariant formu- 
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lation is completely general and applies to all geometries, in- 
cluding irregular shapes, clumpy media, etc. In this general 
form, the path-dependent function ri{y; V) describes the geo- 
metrical morphology, in essence defining the transformation 
from real space to r-space, which might be quite involved. 
In cases of geometrical symmetry, radiative transfer is more 
conveniently formulated with r defined along the axis of 
symmetry. Then ri(y) is defined in an analogous manner, for 
example along the radial direction in the case of spherical 
symmetry (cf Section 



3 HEATED DUST 

The discussion so far was general and applies to all radia- 
tive problems without restrictions. In practice, the scale in- 
variance is mostly useful when the absorption coefficient is 
independent of intensity, and from here on we concentrate 
on such systems. Thus the following discussion excludes line 
emission and photoionization (where the absorption coeffi- 
cient may depend on intensity through its effect on level 
populations) , but is fully applicable to continuum radiation 
- especially by dust, the main thrust of this paper. With 
this restriction, in the absence of a source function radia- 
tive transfer would be described by a homogeneous equa- 
tion. That is, if Ix is a solution when the source term is 
removed then so is alx for any arbitrary constant a. This 
scale invariance is broken by the source term, which makes 
the radiative transfer equation inhomogeneous. In the case 
of a dusty medium, the source function is 



(6) 



where Bx is the Planck function and g{^l',Q) is the angu- 
lar phase function for coherent scattering from direction fi' 
to Cl (e.g. Mihalas 1978). Because of scattering, in general 
Sx contains a dependence on directions. We assumed for 
simplicity single-type dust grains, characterized by a single 
temperature T. The general grain mixture case is similar in 
all fundamental aspects and is discussed separately in ap- 
pendix A. 

Since the scattering term is linear in Ix, it preserves 
the scale invariance of intensity. The only term to break this 
invariance and set a scale for the intensity is the emission 
term, Bx(T). It does that by introducing the dust tempera- 
ture, determined from radiative equilibrium 



dfi / KxSxdX — An 1 KxJxdX 



(7) 



where Jx = J IxdQ/An is the angle-averaged intensity. Both 
sides of this relation vary linearly with the extinction co- 
efficient, therefore the actual magnitude of this quantity is 
irrelevant - only its wavelength variation matters. Indeed, 
Leung (1976) noticed from his numerical calculations that 
T was "rather insensitive" to the absolute values of extinc- 
tion coefficient; our discussion shows that this insensitivity 
is in fact an exact symmetry property of the problem. To 
emphasize this point we write this relation in terms of 



i^x 



(8) 



the extinction coefficient nornmlized to unity at some fidu- 
cial, arbitrary wavelength Aq. Q With the dust source func- 
tion from equation the radiative equilibrium condition 
becomes 



(9) 



q!i.xBx{T)dX= / q^xJxdX 



where gaA = gA(l — t*7A)/(l — roAo) is the absorption ef- 
ficiency normalized to the fiducial wavelength Xq. Because 
the scattering term preserves scale invariance, it does not en- 
ter into this scale-setting equation (which therefore becomes 
meaningless in the case of pure scattering, mx = 1). 

The actual energy input into the dust is external radi- 
ation from a nearby source, such as a star, galactic nucleus, 
etc., of radius rc and intensity I^x- When this source is non- 
spherical (for example disk or torus), the angular profile of 
its geometrical shape, too, must be specified. The dust distri- 
bution has some prescribed arbitrary geometrical shape and 
we denote the distance of its closest point to the source of 
radiation ri. From the previous section, the radiative trans- 
fer problem can not depend separately on rc or ri, only on 
their ratio 



ri 



(10) 



Therefore, the only geometrical property of the radiation 
source relevant to the problem is its angular size at the lo- 
cation of the dust (and its angular shape profile when it is 
non-spherical) . The radiative input to the dust can be char- 
acterized by the external flux at ri , 



F,x{ri) = 



(11) 



This spectral distribution can be specified instead by its 
scale, the bolometric flux at ri 



F,xiri)dX- 



4nrf 



where Lo is the luminosity, and spectral shape 



/c. 



-FcA 



(12) 



(13) 



an r-independent normalized distribution {JfcxdX — 1). 
Then Fcx{ri) ~ /oa-Foi. Consider now the series of mod- 
els in which all properties are held fixed except the input 
bolometric flux Fci. Each value of _Fci uniquely determines 
a corresponding value of Ti , the dust temperature at ri . This 
in turn defines another flux aTi, setting the intrinsic scale 
of dust emission. But the radiative transfer problem can in- 
volve only one radiative scale, therefore the ratio of these 
two fluxes must be a characteristic dimensionless function 
of the model. That is, the ratio 



4aTt 
F,i 



(14) 



is uniquely determined by the dimensionless properties of 
the model such as optical depths, angular shapes and di- 
mensionless density proflle (the factor 4 is introduced for 



Although the absolute value of kx docs not enter independently 
into the radiative transfer problem, it is needed to relate the op- 
tical depth to the total dust mass. 
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convenience). It may also depend on ratios of opacity spec- 
tral averages involving the profile f^x of the input radiation 
and 



(15) 



the normalized (J h\dX = 1) spectral shape of the Planck 
distribution at temperature T\. However, it can not depend 
on any dimensional quantity such as luminosity, linear size 
or density. In Section ^ we present results for in spherical 
symmetry (figure 1) and derive an analytic approximation 
for its r-dependence (equation 

This discussion shows that the luminosity never enters 
into the formulation of the radiative transfer problem, only 
the input flux matters. And rather than specifying the flux 
scale, the dust temperature can be specified instead so that 
flux becomes a derived quantity. A dusty region with pre- 
scribed properties can transport only a single, well defined ra- 
diative flux for a given dust temperature on its heated bound- 
ary. It is determined from the full solution through the scal- 
ing function ^I", which can be viewed as an eigenvalue of the 
model. The fact that the radiative transfer problem can be 
fully specified in terms of temperatures rather than radiative 
flux was first utilized in RR for spherical shells. Equation 
is the general form of this property for all geometries. 

Dust sublimation occurs at some temperature Tsub, de- 
termined by the grain properties. As long as the input 
flux Fei is sufficiently low that T\ < Tsub, the dust tem- 
perature at the inner boundary varies according to eq. |l^ 
(Ti oc _Fei^^*) and different luminosities produce dust emis- 
sion that differs both in overall luminosity and spectral 
shape. This is generally the situation in planetary nebulae. 
However, once _Foi increases to the point that Ti = Tsub, 
the dust temperature cannot rise any more. Instead, further 
increases in luminosity move out the inner boundary of the 
dust so that Fe\ remains fixed at the value that produces 
T\ — Tsub- Thus the temperature profile and emission of the 
dust are set by its internal properties in this regime. Stated 
differently - since F^i uniquely determines Ti, the converse 
is also true and Ti uniquely determines a corresponding Fci . 
Models in which the dust is as close to the radiation source 
as possible have the same Ti and therefore adjust their inner 
boundary to accommodate the external luminosity so that 
they also have the same Fd. Since the bolometric flux is 
now self-regulated by dust sublimation, when the dust is as 
close to the radiation source as possible, the only relevant 
property of the input radiation is its spectral shape fax . 

In this discussion we assumed single-component dust 
and heating by a single external radiative source. Multi- 
component dust mixtures are discussed in Appendix A 
which shows that, although the technical complications in- 
crease, all the essential scaling properties are preserved. The 
same applies to heating by more than one source. Because 
radiative transfer can involve only one intensity scale, the 
addition of an external source can only introduce the ratio 
of corresponding fluxes F^x at every point. 



4 SPHERICAL SYMMETRY 

To further explore concrete consequences of general scal- 
ing we must specify the geometry. Here we take the case of 



spherical symmetry. The essence of most results is preserved 
in other geometries. 

In the case of spherical symmetry, y can be taken as 
dimensionless distance from the center of symmetry, and rj 
and ri{y) defined along the radial path. The radiation source 
is surrounded by a spherical dust shell with inner radius ri 
corresponding to dust sublimation, which we choose as the 
scaling length in equation ^. With this choice, dust sublima- 
tion always occurs aty = l and the dust temperature scale is 
set by T{y=l) = Ti= Tsuh- The dust optical properties are 
fully specified by the wavelength variation of the absorption 
efflciency gaA and the albedo zux- Its spatial distribution is 
described by 



v{y) = 



(16) 



the dimensionless variation of extinction coefficients in the 
radial direction. Note again that ri{y) does not introduce a 
new scale since it is normalized via J^°° ri{y)dy — 1. When 
the dust properties do not vary with distance beyond the 
sublimation point, rj describes the radial law of density vari- 
ation, e.g., y~^, exp{—y), etc. There is no need to specify the 
actual density itself. The extinction scale is set by ro, the to- 
tal optical depth along radial rays at the fiducial wavelength 
Ao. At all other wavelengths, rj — Toqx- At any point y, the 
optical depth from the closest approach to the center along 
a path that makes an angle 9 with the radius vector is 



rx{y,0)^T^ 



y cos 

V 



y2 _J_ yl gjjj^ 



du. 



(17) 



Rowan-Robinson was the first to incorporate geometrical 
scaling for spherical shells into a formalism that specifies 
only scaled radii, and r). 

From the formal solution of the radiative transfer equa- 
tion, the intensity at radius y of a ray inclined at angle 6 to 
the radial direction is 



ix{y,e) = hxe'^^'^'>~>Q[ +iix{y,e). 



(18) 



Here O is the step function (unity for positive arguments, 
zero otherwise) and 



i,x{y,e) = Sxiy',e)e-'^^y'''^-^^^y''^drx{y',e) 



(19) 



is the intensity of the diffuse radiation. By angular integra- 
tion, the angle-averaged intensity is similarly 

Jxiy) = + JM, (20) 

where J^x denotes the diffuse component. We now show that 
Jdx , which is comprised of dust emission and scattering, can 
be expressed explicitly in terms of Jx ■ For the emission term 
this result is straightforward since this term is controlled by 
the dust temperature T, determined from 

qaP(T)T =q^p(Ti)Ti 



J qaxJx{l)d\ 



(21) 



(cf equation Here we introduced the Planck mean of the 
absorption efficiency 



gap(r) = / q^xbx{T)d\. 



(22) 
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For the scattering component, consider first the case of 
isotropic scattering, g(i}',Q) — 1. Then the scattering term 
in the source function is simply ruxJx (eq. and the entire 
diffuse component Jd\ becomes an explicit function of J\. 
Similarly, for the external radiation, radiative equilibrium 
at y = 1 (where T = T\) determines the input bolometric 
flux _Foi, or equivalently ^ (see eq. |l^), as 

(23) 



gap(ri) l-£ 



where 



gac 



]q^xBx{Ti)d\ 



(24) 



Since Jd\ is an explicit function of J\, the same holds for 
5'. Therefore, the entire right-hand-side of equation ^ is ex- 
pressed explicitly in terms of J\ , so this is a self-consistency 
equation for the unknown J\ whose only inputs are dust 
properties and the spectral shape of incoming radiation. As 
an inhomogeneous integral equation, it has a unique physi- 
cal solution. This completes the formulation of the radiative 
transfer problem in accordance with the general scaling dis- 
cussion of the previous section. The only property of the 
external radiation that needs to be specified as input is its 
spectral shape fc\. 

Anisotropic scattering can be handled similarly. In that 
case we write Ix{9) = ^{9)J\, where $ is the intensity an- 
gular distribution, normalized as J $df2 = in. Equation ^ 
for Jx remains the same, only the scattering contribution 
to JdA includes now an integral over <1>. The only technical 
complication is the need to additionally solve the coupled 
equation ( p^ ) for the angular distribution $. The rest of the 
formulation is the same as for isotropic scattering. In com- 
plete analogy, arbitrary geometries can be handled through 
equally straightforward extensions. 

Solution of equation ^ for Jx determines the radial 
variation of the source function in the shell. Once this func- 
tion is known, the flux and surface brightness, the quantities 
relevant for observations, are readily calculated. In particu- 
lar, the overall flux is obtained from the first-moment angu- 
lar integration of equation |l^ which yields 

Fx{y) = %/cAe-""^(^' + Fdxiy), (25) 

similar to the expression for Jx ■ As with the external radia- 
tion, it is convenient to remove the overall scale, contained 
in the bolometric flux F — f FxdX. The spectral shape of 
the overall flux is denoted fx, so that 



Fx{y) = fx(y)F{y), 



F{y) = 4- 
y 



(26) 



Here F\ is the bolometric flux sX y = 1, including both ex- 
ternal and dust contributions, and the second equality ex- 
presses flux conservation. 

The luminosity Le has never entered, and is irrelevant. 
It is important to note, however, that the bolometric flux 
F\ is fully determined even though the overall luminosity 
is not. This flux combines the dimensional scales and 
ri which otherwise do not enter individually, thus dropping 
out of the problem. The corresponding relation, obtained by 
combining equations ^ and can be written as 



where the last equality is written in terms of the source's ef- 
fective temperature Teff- This is the exact expression of a re- 
sult noted in various forms by RR, Laor & Draine (1993) and 
Ivezic & Elitzur (1996a; IE96 hereafter). Therefore, when the 
overall luminosity is known, it can be used to determine the 
radius of the dust sublimation zone. 

The only properties of the central source that affect the 
observed radiation are its spectral shape f^x and angular size 
Sei . The latter enters explicitly only in the expression for the 
intensity, eq. ^ thus it can be determined directly only in 
high-resolution observations that delineate the external con- 
tribution as a bright spot. In most cases, though, the central 
source is unresolved and the quantity relevant for observa- 
tions is the flux. The angular integration that transforms 
intensity to flux removes the explicit dependence on 6^i, as 
is evident from eq. |2^, but the flux may still depend on 6e.i 
because of occultation of the diffuse radiation by the central 
source. Without occultation, the diffuse flux would vanish 
at ri because the intensities of the two streams along any 
ray are then the same (Milne, 1921). Thus the signiflcance 
of occultation can be gauged from the ratio Fd{r\) / F^i of 
the diffuse and external contributions to the bolometric flux 
at ri. Only rays occulted by the central source contribute 
to -FdA(''i), and Bx{Ti) is an upper limit to their intensity. 
Therefore, Fdx{ri) < n6'^iBx{Ti) at all wavelengths and 



^d(ri) 



F,- 



< 



(28) 



Occultation can be ignored whenever Ti <^ Tcs- Comparison 
with eq. ^ shows that this condition is equivalent to 6^1'^ ^ 
1. When this condition is obeyed. 



Fi =Fci, 



fx{y) = fcxe 



-rx(y) 



+ y Fdx{y)/Fi. 



(29) 



and the problem is independent of 9ei. Therefore, as long 
as Tcft 2> Tsub ~ Ti, the only relevant property of the cen- 
tral source is its spectral shape fox', Toff is irrelevant, and in 
the remainder of the discussion we assume that this condi- 
tion is fulfilled JIl When Toff decreases and approaches Tsub, 
a dependence on angular size Oei enters, but Tofi itself re- 
mains irrelevant as long as Ti = Tsub is maintained. Finally, 
scaling breaks down when Tch is so small that Ti < Tsub- 
In this case, T^s must be specified, too, since it determines 
Ti . This breakdown of scaling occurs also when the central 
cavity is dust-free because of factors other then dust subli- 
mation, such as dynamics for example, and Ti < Tsub. 



5 IMPLICATIONS 

When the dust is optically thin at all wavelengths, the solu- 
tion is straightforward. Inserting tx ^ and neglecting the 
diffuse emission, equation ^ gives Jx{y) = Jx{i)/y^ and so 



gap(r) / Tf _ 1 



gap(ri) VTi 



gac 



(30) 



gap(Ti) 

Most of the contribution to the Planck mean comes from 
the distribution peak around 



t For a more detailed discussion of backwarming effects see 
Rowan-Robinson (1982). 
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Xp{T) =4^m 



1000 K 

T ' 



(31) 



so gap(r) ~ q^x (T)- If the corresponding wavelength is 
in the regime where the absorption efficiency drops as a 
power law, gaA oc A"'', we recover the familiar result for the 
dust temperature T = Tiy~'^^^*'^^\ Also, when the spectral 
shape of the external radiation is that of a black-body with 
temperature Te, the optically thin limit of 'i' becomes 



*(r = 0) 



gaAp(To) ^ 
'?aAp{Ti) ' 



(32) 



that is, it is determined by the drop in absorption efficiency 
between the wavelengths corresponding to the temperatures 
of the external source and dust sublimation. 

When the dust is optically thick the problem must be 
solved numerically, and we have developed a computer code 
(DUSTY; Ivezic, Nenkova & Elitzur 1997) to solve the spher- 
ical problem. A description of the numerical procedure is 
provided in Appendix C. Our program solves equation ^ 
for J\{y), with the dust temperature determined from equa- 
tion ^ and the flux scale Fi — AaTi"^ from equation |23[ 
Input includes the optical properties q^x and zux, the dust 
sublimation temperature T^uh , its density distribution ry and 
the external spectral shape f^x- A given set of these quan- 
tities defines a family of solutions, distinguished from each 
other by overall optical depth. Matching a model to obser- 
vations provides the optical depth of the source if there are 
independent estimates for the other properties. Note that 
the only effects on the spectral energy distribution that can 
be meaningfully discerned in observations involve changes of 
spectral features or of slope. 

We present now the results of exact calculations, first 
for what we will use as a "standard" model with = y~^, 
single-component dust grains with size a — 0.05 [im and sub- 
limation temperature T^ub = 700 K, and external radiation 
with the spectral shape of a black body with Tc — 2500 K. 
We then explore the effects of different input quantities on 
the emergent flux and discuss the behavior of the solutions 
with the aid of analytic approximations developed in Ap- 
pendix B. In all cases we vary ry, the overall optical depth 
at 0.55 \im, all the way to 1000 and present solutions for 
both amorphous carbon and silicate grains. Optical proper- 
ties for the former are taken from Manner (1988), for the lat- 
ter from Ossenkopf, Henning & Mathis (1992). The scatter- 
ing is assumed isotropic. Similar explorations of parameter 
space were performed by Leung and by Rowan- Robinson. In 
addition to employing updated optical properties, the cov- 
erage presented here is more extensive and the systematics 
of the presentation more closely reflect scaling. 

Figure 1 displays the variation of ^' with ry. For each 
chemical composition, the optically thin value is a mea- 
sure of the drop in absorption between Ap(ro) ~ 2yim and 
Ap(rsub) — 6 vim (equation ^2|). The thin lines are the an- 
alytic results of equation B7, seen to provide excellent ap- 
proximations to the exact calculations. These analytic ex- 
pressions are not too transparent as they involve spectral 
integrations of the detailed extinction efficiencies. Therefore, 
we discuss the results instead in terms of the gray-opacity 
approximation, equation B9, which maintains most of the 
main features. In particular, it explains the rise of 'I' with 
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Figure 1. Variation of the scaling function ^ (equation 15) with 
overall visual optical depth tv in the "standard model" (see text) . 
Thick lines are the results of exact numerical calculations, the thin 
lines are the analytical approximations of equation B7. 



TV at large optical depths.cJ Indeed, the numerical results 
are well described by the simple analytical fit 



1 -I- 0.005 Tv 



(33) 



^(ry) 
*(0) 

where the power m is 1 for amorphous carbon and 1.25 for 
silicate grains. This provides an adequate approximation for 
^ of the "standard model" over the entire displayed range, 
the accuracy is better than 20% for amorphous carbon, 40% 
for silicates. 

Figure 2 shows the spatial temperature profile for var- 
ious TV. With the density profile i] = y~^, the gray-opacity 
result gives a temperature drop T oc in the in- 

ner regions y ^ tv, followed by a more moderate decline 
T (X y~^^^ at larger radii. Qualitatively, this is the behav- 
ior displayed in the figure, the actual absorption efficiencies 
only modify the details. 

Figure 3 displays the variation of the SED with overall 
optical depth. The shape of the SED is governed by some 
general characteristics, common to all solutions. Because the 
dust cannot be warmer than T^ub, dust emission is negligible 
for A ^ Ap(rsub) and scattered light dominates the diffuse 
component at these wavelengths. As ry increases beyond ~ 
10, the external radiation and scattered light are completely 
absorbed and dust emission dominates at all relevant wave- 
lengths. At each wavelength A, dust radiation is generated 
at radius y where T{y) ~ (40 p.m/A)xlOO K (equation ^l[ ). 
For optically thick configurations, Kwan & Scoville (1976) 



t Formally, this shows that 'I' — > oo when ry — » oo with all other 
properties held fixed, and the inner radius of the shell increases 
without bound (equation ji^ . However, the solutions are derived 
assuming steady-state and the time to reach this stage diverges 
too in that limit. 
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Figure 2. The radial variation of dust temperature in the "stan- 
dard model" for various optical thicknesses, as indicated. 

noted that the SED peaks at the wavelength corresponding 
to Tx ~ 1. At longer wavelengths, the source is always op- 
tically thin and the spectral shape is independent of overall 
optical depth. In this region, the shape of the SED reflects 
only the long-wavelength behavior of the absorption coeffi- 
cient and the density distribution. The net result of all these 
effects is that in optically thick shells, the SED roughly re- 
tains its shape, shifting as a whole to longer wavelengths eis 
TV increases. 

This qualitative discussion adequately explains the gen- 
eral behavior of spectral shapes. We proceed now to explore 

the effects of different input properties on the SED, varying 
them one by one over the entire plausible range; the prop- 
erties not varied are those of the "standard" model. 

5.1 External spectral shape fc\ 

When the external radiation has a black-body spectral shape 
with a temperature Tg, f^x = bx{Tc). figure 4 displays the 
effect of varying Te from 2500 K to 20000 K at various op- 
tical thicknesses. In all cases, longward of 10 ^m Te has 
virtually no effect on the observed spectral shape. The only 
minor effect of Tc, with no practical observational implica^ 
tions, is to slightly modify the relative scale of emission. The 
weak dependence on To was also a result of RR. The reason 
fex h£is no effect at long wavelengths is that the external 




0.1 1 10 100 1000 

X (^m) 

Figure 3. The spectral energy distribution of the "standard 
model" for various optical thicknesses, as indicated. 



contribution is negligible there in comparison with the dust 
emission; although the source is brighter, the dust emission 
covers a much larger surface area. The external radiation 

emerges at short wavelengths because dust emission disap- 
pears at A ^ Ap(rsub). However, even in that spectral region 
the external radiation is visible only when ry 1. At larger 
optical depths it is fully absorbed, becoming irrelevant at all 
wavelengths. Therefore, the external spectral shape can be 
determined only in short-wavelength observations when the 
optical thickness is not too large. This conclusion applies 
to all external spectral shapes, not just those that follow 
Planck's law. 



5.2 Dust sublimation temperature T^ub 

Figure 5 displays the effect of Tsub on the observed spec- 
tral shape. As expected, Taub variations have no effect on 
the shape of the SED at wavelengths shorter than Ap(Tgub)- 
At longer wavelengths, varying Tsub affects only the scale of 
the SED not its shape, similar to the effect caused at those 
wavelengths by variations of Tg. Furthermore, at large opti- 
cal depths (tv ^ 100), spectra produced by varying Tsub are 
indistinguishable at all wavelengths from those produced by 
varying ry. As a result, uncertainty in Tgub translates into 
an uncertainty in the optical depth obtained from a fit to 
the observed SED. 
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Figure 4. Effect of the external radiation spectral shape on the 
emerging SED for various optical thicknesses, as marked. The 
external spectral shape is that of a black-body with temperature 
Te as indicated in the top right panel. In this and all subsequent 
figures, the fixed properties are those of the "standard model" . 

Because the impact of Tsub on the spectral shape is mi- 
nor and poorly separated from other effects, Tsub cannot be 
well determined from the SED. The only regime where such 
determination can be realistically attempted is intermediate 
optical depths, tv ^ 10. In general, the dust sublimation 
temperature is more reliably determined in high-resolution 
observations whenever they are capable of resolving the dust 
condensation zone. From equation 

^Tt = (34) 

where -Fobs and Si^obs are, respectively, the observed flux and 
dust condensation zone angular diameter (see also IE96). 
Whenever feasible, this is the more practical method and 
we have successfully employed it for analysis of IRC-l-10216 
(Ivezic & Elitzur, 1996b). 

5.3 Dust distribution 

Figures 6-8 display the effect of the dust density distribu- 
tion, presenting the results for power-law variation r] oc y^^ 
for various non-negative p. As is evident from these flgures, 
for given dust optical properties r) is the input quantity with 
the most substantial impact on the long wavelength emis- 
sion, the only one to produce an unambiguous observational 
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Figure 5. Effect of dust sublimation temperature Tgub, indicated 
in the top right panel, on the emerging SED. 

signature. Different density distributions can be easily dis- 
tinguished by colors involving A J> 10 \im. 

Varying the dust density distribution affects the emerg- 
ing spectrum in two ways. It changes the amount of dust re- 
sponsible for the emission at a given wavelength and it also 
modifies the temperature profile. However, the differences 
between the temperature profiles play only a secondary role 
to those in the matter distribution itself. The flux at wave- 
length A is dominated by emission from radius y where the 
dust temperature T obeys Ap(r) ~ A, so redistributing the 
material between the various radii directly modifies the spec- 
tral shape. Power-law distributions with p ^ 1 must include 
a cutoff to ensure finite column density. We therefore make 
a distinction between two families of distributions, discussed 
separately. 

5.3.1 Steep distributions; p> 1 

Steep density distributions do not require a cutoff. Fig- 
ure 6 displays the SEDs for p = 2, 3/2 and 1 + e, where 
e <C 1. The spectral shapes become flatter at long wave- 
lengths as p decreases since for the same overall optical 
depth, relatively more material is placed at larger radii (i.e., 
low temperatures). In particular, for power-law variations 
gaA oc A"'' and r) oc ^ Harvey et al. (1991) show that 
A/a oc \-(0+'^)iv-'^)/'2 in the optically thin regime. We find 
this relation useful for qualitative understanding of the SED 
past its peak. Since the differences among SEDs reflect dust 
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Figure 6. Effect of dust density profile r/ on the emerging SED 
for power-law profiles steeper than , indicated in the top right 
panel. 

emission, these differences disappear at A ^ A(Tsub)- Indeed, 
as the figure shows, the spectral shapes are indistinguishable 

at such wavelengths, therefore the density distribution can- 
not be determined from the observed SED at A ^ A(rsub). 
At large optical depth, rv ^ 50, the emission peak shifts to 
wavelengths longer than A(Tsub) and the SEDs differ from 
each other on both sides of the peak. 

5.3.2 Flat distributions; p ^ 1 

Power-law distributions with p ^ 1 require a cutoff j/out 
to ensure a finite column density, introducing an additional 
independent parameter, the shell relative thickness Y = 
Uout — 1- This introduces in turn a lower limit for the temper- 
ature and a corresponding wavelength Aout = Ap(T(yout)). 
Dust emission is mostly confined to wavelengths shorter that 
Aout, which thus serves as an effective long- wavelength cutoff 
for the SED. 

The cutoff effect is best seen by comparing the SEDs 
for ?7 oc of figure 7 to those displayed in figure 6 for 
77 oc j/^^^^'-*. The latter density distribution does not require 
a y-cutoff and displays the intrinsic effect of this power law 
for an infinite shell. The strict distribution must include 
a cutoff, which is varied in figure 7 over a wide range. For all 
practical purposes, the SED for Y = 10^ is indistinguishable 
from that for the infinite shell. The difference between the 
two is noticeable only at wavelengths longer than ^ 300 
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Figure 7. Effect of the shell's geometrical thickness on the emerg- 
ing SED for dust density profile r) (xy~^ . Only the relative thick- 
ness Y = Ar/ri, indicated in the top right panel, is relevant 

[im, where the SED of the finite shell dips down. As the 
cutoff Y decreases, Aout decreases too as more cold dust is 
removed from the shell and the decline of the SED is pushed 
to shorter wavelengths. 

Figure 8 presents the SEDs for constant-density shells 
with various thicknesses. The gray-opacity result of equation 
B9 yields 

Aout oc [STviY + lj + iY + lfY^" (35) 

for such shells when they are not too thin either physically 

(Y ^ 1) or optically (rv ^ 1). This simple relation explains 
adequately the behavior of the long-wavelength cutoff evi- 
dent in the various panels of figure 8. When Y <^ ry, the 
cutoff wavelength varies in proportion to {ryY)^^''^. When 
Y ^ Tv, the cutoff wavelength becomes independent of 
rv and varies in proportion to Y^^^. In geometrically-thin 
(Y <C 1) optically-thick (ry S> 1) shells, the cutoff wave- 
length varies in proportion to tv^^*, independent of the shell 
geometrical thickness. 

It is important to note that, although steep density dis- 
tributions do not require a formal cutoff, in practice such 
cutoffs always exist because of the finite sizes of dust shells. 
These finite sizes affect the long wavelength part of the SED 
in exact analogy to the cutoff of a flat distribution. The fi- 
nite size of a steep density distribution is relevant only for 
observations at wavelengths longer than the corresponding 
Aout- In all cases, the cutoff effect can be easily masked by 
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Figure 8. Same as figure 7, but for a constant density distribu- 
tion. 



Figure 10. Effect of grain size, indicated in the figure, on the 
emerging SED. 
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Figure 9. Variation of the absorption efficiency spectral shape 
with grain radius, indicated in the right panel. 

background emission because it generally involves low tem- 
peratures. 

5.4 Optical properties; grain chemistry and size 

The optical properties gaA and vjx depend on the grain 
chemical composition and size. The effect of composition 
on the SED can be discerned in all figures by comparing 
the panels for amorphous carbon and silicate grains. The 
significance of grain size a can first be gauged from figure 
9. For the two chemical compositions and various radii a 
it displays QaA , obtained from a Mie theory calculation with 



the appropriate dielectric constants. Except for spectral fea^ 
tures, the extinction efficiency is a self-similar function of 
the scaled variable 2na/\, approximately constant at short 
wavelengths and switching to a power-law decline with in- 
creasing A at long wavelengths. In this approximation, the 
shape of <7aA is independent of grain size for a ^ A/27r. There- 
fore, at wavelengths longer than visual, grain sizes are irrel- 
evant as long as a ^ 0.05 \iTa. On the other hand, larger 
grain radii will significantly affect the results because they 
are effectively equivalent to selecting grains with different 
radiative properties. Increasing the grain size beyond 0.05 
^i.m amounts to increasing gaA at every wavelength longer 
than visual. Therefore, for the same rv, increasing the grain 
size increases all far-IR optical depths. 

Figure 10 displays the SEDs for our "standard model" 
with grains of different sizes and various optical depths. The 
results for all grains smaller than 0.05 \un should be the 
same as for a = 0.05 [im. As could be expected, the effect 
of grain size becomes significant only for a > 0.1 |i.m. The 
difference between the SEDs for a = 0.05 and 0.1 (im is gen- 
erally small, especially for amorphous carbon. Larger grains, 
a ^ l|-im, generate some structure at intermediate optical 
depths. This structure, evident in the top two panels for 
amorphous carbon, arises purely from optical depth effects, 
reflecting the interplay between the stellar component and 
dust emission and scattering. It can give a false impression 
of dust spectral features, although they do not exist in the 
actual absorption efficiency (see figure 9). 
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6 DISCUSSION 

Scaling has far reaching consequences for modeling and anal- 
ysis of IR emission from radiatively heated dust. Thanks to 
scaling, the number of independent input properties is signif- 
icantly reduced. Only two scales need be specified for a com- 
plete solution — the overall optical depth ry and the dust 
sublimation temperature Tsub. All other input properties in- 
volve only dimensionless, normalized profiles. For the exter- 
nal radiation, only the spectral profile /ca is required. For 
the dust density, only the normalized distribution rj (equa- 
tion ^ is needed. And for the dust optical properties, only 
the spectral shapes of absorption and scattering efficiencies 
enter. 

Scaling enables a meaningful systematic study of the 
observational impact of each independent input property on 
the emerging IR radiation. The results of the systematic cov- 
erage of parameter space presented in the previous section 
outline the range of spectral shapes produced by spherical 
configurations. They show which properties delineated in 
observations can and cannot be explained in this symmetry. 
For example, the IR spectral shape XFx oc A^^'''' is some- 
times used as an indicator of disk geometry (e.g. Hillenbrand 
et al. 1992). In fact, this spectral behavior is also produced 
by spherical distributions with rj oc with p ~ 1.5 if they 
are optically thin in the wavelength region where qx oc A"*^, 
since j3 is typically 1-2 (Miroshnichenko, Ivezic & Elitzur 
1997). 

Our results suggest a natural classification scheme for 
IR spectra. For given grain optical properties, every rj pro- 
duces a distinct family of solutions, with position within the 
family determined by rv. For the most part, T^uh and f^x 
have only a minor effect. Therefore, the SEDs of all astro- 
nomical objects that share the same dust density distribu- 
tion become a one-parameter family in which each member 
is fully characterized by its tv- This leads to correlations 
among spectral properties and structure in color-color dia- 
grams. Scaling breaks down when the dust inner boundary 
is controlled by effects other than dust sublimation. In such 
sources, the inner-boundary dust temperature Ti (< Tsub) 
becomes an additional independent parameter, varying with 
the external fiux according to Ti oc F^i^''^. 

Various classes of astronomical objects show evidence 
for scaling, indicating that in these sources the dust in- 
ner boundary is controlled by sublimation. A common spec- 
trum for compact HII regions was noted in Rowan- Robinson 
(1979) and modelled in RR. This was further emphasized 
by Wood and Churchwell (1989), who point out that for all 
ultra-compact HII regions observed by IRAS, "the shapes 
of flux density distributions are strikingly similar", a fact 
also recognized by Chini et al. (1986). The trends observed 
among spectral properties of IR radiation from young stellar 
objects elicited a number of classification schemes. Adams, 
Lada & Shu (1987) employed the slope of the SED, Adams 
(1990) the visual extinction derived from model spectra, 
Ladd et al. (1991) the mean frequency of the spectrum and 
Myers & Ladd (1993) the corresponding bolometric tem- 
perature. Myers & Ladd also noted that all these schemes 
are equivalent, but the reasons for the trends and for this 
equivalence were not clear. Scaling provides a simple ex- 
planation for both. Furthermore, various distributions pre- 
sented in Myers & Ladd show that IR spectral properties 



are indeed independent of luminosity, as prescribed by scal- 
ing. Finally, the structure in color-color diagrams implied by 
scaling had been noticed in the IRAS colors of a variety of 
objects. These include T Tau stars (Harris, Clegg & Hughes, 
1989), BN objects (Henning, Pfau & Altenhoff, 1990), HII 
regions (Hughes & MacLeod, 1989; Wood & Churchwell, 
1989), H2O masers in star- forming regions (Wouterloot & 
Walmsley, 1986) and late-type stars (van der Veen & Habing, 
1988). Various scaling properties of the latter were noted 
previously by Rowan- Robinson and Harris (1983 and refer- 
ences therein). 

Dusty winds around late-type stars provide an even 
tighter form of scaling (IE95, IE96). The structure of these 
winds is controlled by radiation pressure on the dust grains, 
therefore the density distribution rj need not be prescribed 
in advance. Instead it is determined from the hydrodynamics 
coupled to the radiative transfer, and the solution does not 
require any other parameter in addition to optical depth. 
Therefore, IR emission from all late-type stars with a given 
dust composition can be fully classified in terms of a single 
parameter, rv- Detailed analysis shows that the data fully 
corroborate all the correlations predicted by scaling for these 
sources. We are in the midst of a study that extends this de- 
tailed analysis to the IRAS data for all the Galactic point 
sources. The scaling correlations displayed by these data will 
enable us to identify the dust density distributions of vari- 
ous classes of Galactic objects. The results of this study will 
be reported separately in a forthcoming publication. 
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A convenient indexing of the dust components is by 
decreasing sublimation temperature. Consider a spherically 
symmetric configuration with a prescribed radial variation of 
overall dust density when the shell inner radius is controlled 
by sublimation of the first component at temperature Tsub,!. 
The problem is identical to that of a single-component dust 
as long as the temperature of component 2 at dimensionless 
radius y, determined from equation h.4, remains higher than 



its sublimation temperature rsub,2. The radius j/2 at which 
component 2 is included in the mix is thus determined from 
radiative transfer, and so on. The only difference in the solu- 
tion procedure from the case of single-type dust is that the 
dimensionless density profile rj[y) cannot be fully prescribed 
beforehand since the determination of the sublimation radii 
1/2, 1/3, etc., requires solution of the radiative transfer prob- 
lem. Scaling is preserved but the exact form of r](y) requires 
an iterative procedure. The extension of this procedure to 
arbitrary geometries is straightforward. 



APPENDIX A: GRAIN MIXTURES 

When the dust is composed of different types of grains, its 
extinction coefficient is 

n\ = Pd}Pi — -Q\i- (Al) 
^ nii 

i 

Here pd is the overall dust density and Pi is the fraction of 
that density due to the i-th component (^Pi = 1). The 
mass, radius and extinction efficiency of a single grain from 
this component are denoted m^, and Q\i, respectively. 
The different components refiect possible variations of chem- 
ical composition and/or grain size. For simplicity we treat 
the components as discrete; generalization to the continuous 
case is trivial. With the aid of equation ^ the dust source 
function is 



4n 



(A2) 



where u^xi and gi are, respectively, the albedo and angular 
phase function of the i-th component and 



Xxi 



Pi{Tvai/mi)Qxi 
'^^Pi{-naymi)Qxi 



(A3) 



is the weight of its contribution to the extinction coefficient 
at wavelength A. Each component may have a different tem- 
perature, Td,i, determined from overall fiux conservation. 
This equilibrium relation becomes a vanishing sum of the 
fractions pi with coefficients that depend on the correspond- 
ing individual temperatures Td,i- Since the temperature of 
a dust component cannot depend on its relative abundance, 
each coefficient in the sum must vanish separately. This pro- 
duces the temperature equation of the i-th component 



(A4) 



(?aA,i-BA(Td,i)rfA = / q^xsJxdX 



where q^x,i = gAi(l - ■s7Ai)/(l - c^Aoi) is the absorption 
efficiency of the i-th component normalized to the fidu- 
cial wavelength Aq, analogous to equation ^ for a single- 
component dust. 



APPENDIX B: APPROXIMATE SOLUTIONS 
IN SPHERICAL SYMMETRY 

The radial variation of temperature requires the spatial pro- 
file of Jx, equivalent to the spatial variation of both the 
overall radiative energy density J[y) = J Jx{y)dX and the 
frequency profile Jx/J- The former can be obtained from 
the wavelength-integrated radiative moment equations when 
they are closed by some suitable approximation. Auer (1984) 
proposed the closure relation 



3P = 4nJ + 2F, 



(Bl) 



where P is the bolometric radiation pressure, obtained from 
the second angular moment of the intensity. This relation 
improves on the classic Eddington approximation (3P = 
47tJ) by providing the correct behavior in both the opti- 
cally thin and thick limits. With this approximation, the 
radiation-pressure equation becomes 



47t : 

dy 
where 



1,2 



- + 3ToqF{y)v{y) 

y 



9F = / <l\f\d\ 



(B2) 



(B3) 



is the extinction efficiency averaged with the fiux spectral 
shape. Solution of this equation requires a boundary con- 
dition. We assume that the energy density vanishes when 
y oo, setting aside for a moment the case of finite shells. 
Then the solution is 



4nJ{y) = Fi 



+ 3ro 



, , , \dy 

qv[y)v[y)-^ 
y 



(B4) 



This result requires the fiux-averaged gp, which cannot be 
calculated without the spectral shape of the diffuse com- 
ponent of the ffux (see equation p9[ ). We will approximate 
this with the Planck spectral shape at the local tempera- 
ture, namely F^x oc bx{T). This can be expected to give an 
adequate approximation at large optical depths, less so at 
small ones, where the significance of the diffuse component 
is diminished anyhow. With this approximation, 



QFiy) = / qxfexe 



'Tx(y) 



d\ 
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+ gp(r(y)) 1- / Ue 



-Tx(y) 



(B5) 



This completes the approximation for J. Next we need the 
spectral shape of the diffuse component J^a- This can be 
reasonably expected to resemble the spectral shape of the 
diffuse flux -FdA, so we approximate it, too, with bx{T). Then 
the radiative equilibrium equation (H) gives 



4nJ{y) = 4aT^ + ^ //oAe"^^'"' 



'7ap(r) 



dX. 



(B6) 



Equations B4 and produce the temperature profile as a 
function of y, with Fi, or equivalently 'I' (see equation 
obtained from the expressions at y = f . The results are 



gap(r) 



+ 3ro / qF{y)'n{y)%, 

J y y 
-^^,+3ro r qHy)v{y)% (B7) 



note that equation |30| is recovered when tq = 0. Together 
with equation B5, this completes the solution for the tem- 
perature profile for any given density profile rj that ex- 
tends to infinity. The solution is only in implicit form be- 
cause T enters on the right-hand-side through the Planck- 
averages qp{T) and qa,p{T). Therefore the complete solution 
requires iterations, but these are straightforward and con- 
verge rapidly. We find that these analytic results provide 
an excellent approximation for the exact solutions, typically 
within 10%. 

When the shell is finite, a boundary condition must be 
specified for its energy density at the outer edge j/out- From 
the definitions of J and F as angular moments of the inten- 
sity it follows that 



7-Fi 

47tJ(yout) = ^— , 

Vout 



(B8) 



where 7 varies from 1 when the angular distribution of the 
emergent radiation is directional to 2 when it is isotropic. 
Leaving 7 as an unknown factor, this term is ju st a dded to 

simply 



B4 



the solution on the right-hand-side of equation 
adding 7/yout to the right-hand-side of the two relations in 
equation B7. So the analytic approximation for finite shells 
includes a term that remains uncertain to within a factor 
of 2. However, this term is significant only in shells that 
are thin both physically and optically (yout — 1, ro < 1) 
and our numerical results show that 7 ~ 1.5 is an adequate 
approximation in virtually all cases. 

A different approximate solution scheme was proposed 
in RR. That scheme does not fully incorporate scattering but 
its iterations include a variable Eddington factor, therefore 
it is in principle capable of convergence to the exact solu- 
tion in the absence of scattering. Since our scheme involves 
a prescribed Eddington factor (eq. Bl) it only converges to 
an approximate solution, not the exact one. However, when 
scattering is important it has the edge and we usually obtain 
a level of accuracy comparable to RR with siginficantly fewer 
computations, mostly because we do not have angular inte- 
grations. Furthermore, our approximate solution becomes 
fully explicit in the "gray opacity" case, the approximation 
that assumes frequency independent extinction. Then qx = 



1 and the solution reverts to the simple, explicit form 
1 



rr-f \ 4 

*x ( — 



* = 1 + 3ro 



y2 



(B9) 



In the case of a finite shell, the right-hand-sides of both 
equations again are augmented by j/yout- This gray-opacity 
result follows directly from the Auer closure relation and 
does not require assumptions about the spectral shape of 
the diffuse component. It provides a closed-form expression 
that readily explains many properties of the full solutions, 
as pointed out in section 5. 



APPENDIX C: NUMERICAL PROCEDURES 

The code DUSTY, which solves the spherical radiative 
transfer problem employing the scaling approach described 
here, is publicly available for general use (Ivezic, Nenkova & 
Elitzur 1997). The numerical calculations are conveniently 
performed with the aid of the dimensionless energy density 



4ny 



Jx 



(CI) 



Extracting the y'^ radial dilution reduces the dynamic range 
of the density, helping the numerical accuracy. For isotropic 
scattering, ma at radius y obeys the integral equation 

/oAe-^^(«) 



ux(y) 



1 

+ 2 



zuxux{y') 



(l-t^A)*( J 



hx{T') 



^^Tx(y',^)-Tx(y,^) drxiy',^I)d^,, (C2) 

where /i = cos 9. In this equation T' — T(y') is obtained 
from 



gaP (T') fT'Y _ J q^xux{y')d\ 
Ti 



gap(ri 

and from 
1 



gaP (Tl 



/ gaAMA(l)dA 



gaAMA(l)rfA. 



(C3) 



(C4) 



Once a radial grid is prescribed, the numerical evaluation 
of the integral of any radial function is transformed into 
multiplication with a matrix of weight factors determined 
purely by the geometry (cf Schmid-Burgk, 1975). Therefore, 
if T{y) and are given, ux can be obtained for all radial 
grid points through matrix inversion. This provides a direct 
solution to the full scattering problem, obviating the need 
to iterate over ux{y) itself. 

We start the calculation with T{y) and ^ from the an- 
alytic solution of appendix B, calculate the corresponding 
ux{y) and iterate until convergence. Because the scattering 
problem is solved directly by the matrix inversion, conver- 
gence is very rapid; an accuracy of 10"** is achieved with 
fewer than 10 iterations even for the largest optical depths 
considered here. The solution is accomplished once the den- 
sity ua does not vary, within a prescribed tolerance, when 
the grid density is increased. The number of radial grid 
points we needed varied from ~ 20 for moderate optical 
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depths (rv ^ 10) to 80 for rv = 1000. The number of an- 
gular grid points used in the ^-integration is typically 2-3 
times the number of radial grid points. Once the energy den- 
sity u\ is found, the flux spectral shape f\{y) is calculated 
from the same integral with d/x replaced by iidfj,. 

We have compared our code with a number of bench- 
mark solutions and verified its accuracy to better than 0.1% 
(Ivezic et al. 1997). Most of the computer time is spent on 
calculation of the weight matrix for spatial integrations. This 
matrix is determined purely by geometry and does not vary 
between iterations, giving our method an edge in calcula- 
tion speed over other schemes (e.g. iterations over variable 
Eddington factor). In particular, large optical depths do not 
increase the computational difficulties as much as they do 
in purely iterative schemes; we encountered no difficulties 
with rv as large as 2000. In addition, the inclusion of scat- 
tering iu the direct matrix inversion eliminates the problems 
that large albedos present in other schemes. Our method re- 
mains stable even for albedos that approach unity at some 
wavelengths. 
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